{ "cells": [ { "cell_type": "markdown", "id": "3906abf2", "metadata": {}, "source": [ "---\n", "format:\n", " html:\n", " other-links:\n", " - text: This notebook\n", " href: L4 Solving nonlinear equations in 1d.ipynb\n", "---" ] }, { "cell_type": "code", "execution_count": 2, "id": "81afe21c", "metadata": {}, "outputs": [], "source": [ "using Plots\n", "using LaTeXStrings" ] }, { "cell_type": "markdown", "id": "aad6bc09-0de2-4531-9043-0bacc8e41409", "metadata": {}, "source": [ "# Solving nonlinear equations in 1d" ] }, { "cell_type": "markdown", "id": "c79162dc", "metadata": {}, "source": [ "
\n", "\n", "Note. This chapter will be mainly based on Chapter 1 of An Introduction to Numerical Analysis by Süli, Mayers. These notes are mainly a record of what we discussed and are not a substitute for attending the lectures and reading books! If anything is unclear/wrong, let me know and I will update the notes.\n", "
" ] }, { "cell_type": "markdown", "id": "bdb3c1d1", "metadata": {}, "source": [ "
\n", "\n", "Note: \n", "These notes contain some (but not all) of the content of the lecture. Please let me know if there are any additions that you feel are particularly important and let me know if you find any errors. \n", "
" ] }, { "cell_type": "markdown", "id": "316c2311", "metadata": {}, "source": [ "## Previously..." ] }, { "cell_type": "markdown", "id": "43fab67d", "metadata": {}, "source": [ "
⚠ Subtractive cancellation \n", "\n", "When adding or subtracting numbers so that the result is much smaller in magnitude. \n", "\n", "Example 1: We know $x = -1.0000000000000000000147$ to $23$ significant digits but $x + 1 = -1.47 \\times 10^{-20}$ to only $3$ significant digits \n", "\n", "Example 2: Quadratic formula when $-b \\approx \\pm \\sqrt{b^2-4ac}$ \n", "\n", "Example 3: $f(x) = \\frac{e^x - 1}{x}$ vs $p(x) = 1 + \\frac{x}{2!} + \\frac{x^2}{3!} + \\dots + \\frac{x^8}{9!}$\n", "\n", "...\n", "
" ] }, { "cell_type": "markdown", "id": "7da07edb", "metadata": {}, "source": [ "Relative condition number:\n", "\n", "\\begin{align}\n", " \\kappa_f(x) &= ... = \\left| \\frac{x f'(x)}{f(x)} \\right|\n", "\\end{align}" ] }, { "cell_type": "markdown", "id": "55eab8b9", "metadata": {}, "source": [ "
Definition. \n", "\n", "If the condition number is large, we say the problem is *ill-conditioned* or the problem is badly conditioned \n", "When the error can not be explained by an ill-conditioned problem, we will say the algorithm is *unstable*.
" ] }, { "cell_type": "markdown", "id": "a46bc6c8", "metadata": {}, "source": [ "## Chapter 2\n", "\n", "* Aim: for $f: [a,b] \\to \\mathbb R$, solve $f(\\xi) = 0$ for $\\xi \\in \\mathbb R$, \n", "* Also, can find $\\xi$ such that $f(\\xi) = c$ (by considering $x \\mapsto f(x) - c$) or stationary points (by considering $f = F'$),\n", "* Change of sign theorem. Converse is not true\n", "* Can rewrite as $g(x) = x$ (in multiple ways)\n", "* Brouwer's fixed point theorem" ] }, { "cell_type": "markdown", "id": "02cd7482", "metadata": {}, "source": [ "
Example. \n", "\n", "Solve $f(x) = x$ where $f(x) := e^x - 2x - 1$. \n", "Apply change of sign theorem to $f$ to conclude there exists $\\xi \\in [1,2]$ such that $f(\\xi) = 0$.\n", "
" ] }, { "cell_type": "code", "execution_count": 3, "id": "76076758", "metadata": {}, "outputs": [], "source": [ "# definitions for the example f(x) = e^x - 2x - 1\n", "\n", "function simple_iteration( g, x1, N=100, tol=1e-10 )\n", " x = [ x1 ]\n", " for n in 2:N\n", " push!( x, g(x[n-1]) )\n", " if (abs(g(x[end]) - x[end]) < tol)\n", " break\n", " end\n", " end\n", " return x\n", "end\n", "\n", "f = x -> exp(x) - 2*x - 1;\n", "g = x -> log( 2*x + 1 ); \n", "h = x -> (exp(x) - 1)/2;\n", "id = x -> x;\n", "\n", "f_prime = x -> exp(x) - 2;\n", "f_2prime = x -> exp(2);\n", "g_prime = x -> 2/( 2*x + 1);\n", "\n", "ξ = simple_iteration(g, 1.0, 100, 1e-15)[end];\n", "ζ = 0.0;" ] }, { "cell_type": "code", "execution_count": 4, "id": "5f8b44f0", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(f(1), f(2)) = (-0.2817181715409549, 2.3890560989306504)\n" ] }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@show ( f(1), f(2) )\n", "plot( f, 1, 2, label=L\"y = f(x)\" )\n", "hline!([0], label=L\"y = 0\", linestyle=:dash)\n", "scatter!( [1 2], [f(1) f(2)] , label=[L\"f(1)\" L\"f(2)\"]) " ] }, { "cell_type": "markdown", "id": "e1db43ec", "metadata": {}, "source": [ "
\n", "\n", "We could solve $f(x) = 0$ by instead considering the fixed point problems $g(x) = x$ or $h(x) = x$ where\n", "\n", "\\begin{align}\n", " g(x) = \\log( 2x + 1 ), \\quad \\text{and} \\quad h(x) = \\frac{e^x - 1}{2}\n", "\\end{align}\n", "\n", "Brouwer's fixed point theorem applies to $g$ but **not** $h$ on $[1,2]$:\n", "
" ] }, { "cell_type": "code", "execution_count": 5, "id": "be89f658", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(g(1), g(2)) = (1.0986122886681098, 1.6094379124341003)\n", "(h(1), h(2)) = (0.8591409142295225, 3.194528049465325)\n" ] }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@show (g(1), g(2));\n", "@show (h(1), h(2));\n", "\n", "plot( [g, id, h], 1, 2, \n", " label=[L\"g\" L\"x\" L\"h\"], lw=2,\n", " color=[\"red\" \"black\" \"green\"], linestyle=[:solid :dash :solid] )\n", "hline!( [1,2], linestyle=:dash, primary=false, color=\"black\" )" ] }, { "cell_type": "markdown", "id": "6449418f", "metadata": {}, "source": [ "* Lemma: if $(x_{n}) \\to \\xi$ then $\\xi$ is a fixed point of $g$. Proof: Use the fact that $x_{n+1}\\to \\xi$ and $g$ is continuous so $g(x_n)\\to g(\\xi)$.\n", "\n", "
Example (cont.). \n", "\n", "Returning to $f(x) = e^x - 2x - 1$, do the iterative methods $x_{n+1} = g(x_n)$ and $y_{n+1} = h(y_n)$ converge for $x_1 = y_1 = 1$?\n", "\n", "
" ] }, { "cell_type": "code", "execution_count": 6, "id": "1918945a", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "22-element Vector{Float64}:\n", " 0.001\n", " 0.001998002662673058\n", " 0.0039880425020066826\n", " 0.007944444173298813\n", " 0.0157639813122604\n", " 0.031041161862191746\n", " 0.060231437446657\n", " 0.11374188108461382\n", " 0.2049663522236682\n", " 0.3435419759087765\n", " 0.5230015662957427\n", " 0.7158881986091449\n", " 0.8886220179437329\n", " 1.0214590819888538\n", " 1.1128169773483771\n", " 1.1711295060291513\n", " 1.206646929282068\n", " 1.2276777660025902\n", " 1.2399253555448633\n", " 1.2469893937540422\n", " 1.251041140573826\n", " 1.25335772906099" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = simple_iteration( g, 1e-3, 22 )" ] }, { "cell_type": "code", "execution_count": 7, "id": "f4dbff2e", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[36m\u001b[1m[ \u001b[22m\u001b[39m\u001b[36m\u001b[1mInfo: \u001b[22m\u001b[39mSaved animation to c:\\Users\\jack6\\Math 5485\\Pictures\\simple_iteration_stable.gif\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "Plots.AnimatedGif(\"c:\\\\Users\\\\jack6\\\\Math 5485\\\\Pictures\\\\simple_iteration_stable.gif\")" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot( [id, g] , min(-0.25,x[1]-1e-3), max(2,x[1]+1e-3), \n", " label=[ L\"x\" L\"g(x)\"], title=L\"\\textrm{simple~iteration~} x_{n+1} = g(x_{n})\",\n", " color=[\"black\" \"red\"], linestyle=[:dash :solid], lw=2)\n", "\n", "anim = @animate for n ∈ 2:2(length(x)-1)\n", " if (n%2 == 0)\n", " k = Int(n/2)\n", " plot!( [x[k], x[k]], [x[k], x[k+1]], \n", " primary=false, lw=2, color=\"blue\")\n", " else\n", " k = Int((n+1)/2)\n", " plot!( [x[k-1], x[k]], [x[k],x[k]], \n", " primary=false, lw=2, color=\"blue\")\n", " end\n", "end\n", "gif(anim, \"Pictures/simple_iteration_stable.gif\", fps = 3)" ] }, { "cell_type": "code", "execution_count": 8, "id": "7e42e536", "metadata": {}, "outputs": [], "source": [ "y = simple_iteration( h, ξ-1e-2, 15 ); " ] }, { "cell_type": "code", "execution_count": 9, "id": "294c91b2", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[36m\u001b[1m[ \u001b[22m\u001b[39m\u001b[36m\u001b[1mInfo: \u001b[22m\u001b[39mSaved animation to c:\\Users\\jack6\\Math 5485\\simple_iteration_unstable.gif\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "Plots.AnimatedGif(\"c:\\\\Users\\\\jack6\\\\Math 5485\\\\simple_iteration_unstable.gif\")" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot( [id, h] , min(-0.25,y[1]-1e-3), max(2,y[1]+1e-3), \n", " label=[ L\"x\" L\"h(x)\"], title=L\"\\textrm{simple~iteration~} y_{n+1} = h(y_{n})\",\n", " color=[\"black\" \"green\"], linestyle=[:dash :solid], lw=2)\n", "\n", "anim = @animate for n ∈ 2:2(length(y)-1)\n", " if (n%2 == 0)\n", " k = Int(n/2)\n", " plot!( [y[k], y[k]], [y[k], y[k+1]], \n", " primary=false, lw=2, color=\"blue\")\n", " else\n", " k = Int((n+1)/2)\n", " plot!( [y[k-1], y[k]], [y[k],y[k]], \n", " primary=false, lw=2, color=\"blue\")\n", " end\n", "end\n", "gif(anim, \"simple_iteration_unstable.gif\", fps = 3)" ] }, { "cell_type": "code", "execution_count": 10, "id": "1ea03f4f", "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot( [g, h] , -0.25, 2, \n", " label=[ L\"g(x)\" L\"h(x)\"], color=[\"red\" \"green\"], lw=2)\n", "plot!( x, g.(x), seriestype=:scatter, label=L\"x_{n+1} = g(x_n)\" ,color=\"red\")\n", "plot!( y, h.(y), seriestype=:scatter, label=L\"y_{n+1} = h(y_n)\", color=\"green\" )" ] }, { "cell_type": "markdown", "id": "b047b5b1", "metadata": {}, "source": [ "
\n", "\n", "Here, we see that $x_n \\to \\xi \\approx 1.26 \\in [1,2]$ but $y_{n} \\to 0$ so the formulations $g(x) = x$ and $h(x) = x$ are equivalent, but they lead to different numerical methods.\n", "
" ] }, { "cell_type": "markdown", "id": "3a0e4661", "metadata": {}, "source": [ "
Definition. \n", "\n", "Let $(x_n) \\to \\xi$. Suppose that there exist $\\alpha\\geq1$ and $\\mu \\geq 0$ such that\n", "\n", "\\begin{align}\n", " \\lim_{n\\to\\infty} \\frac{\\left|x_{n+1} - \\xi\\right|}{\\left|x_n - \\xi\\right|^\\alpha} = \\mu\n", "\\end{align}\n", "\n", "* If $\\alpha = 1$ and $\\mu \\in (0,1)$, we say the convergence is *linear*,\n", "* If $\\alpha = 1$ but $\\mu = 0$, we say the convergence is *superlinear* (i.e. faster than linear),\n", "* If $\\alpha = 1$ and $\\mu = 1$, we say the convergence is *sublinear* (i.e slower than linear),\n", "* If $\\alpha = 2$, we say the convergence is *quadratic*,...\n", "\n", "In the case where *(i)* $\\alpha = 1, \\mu \\in (0,1)$ or *(ii)* $\\alpha>1, \\mu> 0$, we say\n", "* $(x_n) \\to \\xi$ with *order* $\\alpha$, \n", "* $\\mu$ is the *asymptotic error constant*, and\n", "* $\\rho := - \\log_{10} \\mu$ is the *asymptotic rate of convergence*.\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "8ba2e2cc", "metadata": {}, "source": [ "* Define the errors $e_n \\coloneqq |x_{n}-\\xi|$. Then, we have $e_{n+1} \\approx \\mu e_{n}^\\alpha$,\n", "* For linear convergence, $e_{n+N} \\approx \\mu^n e_N = 10^{-\\rho n} e_N$ and so $\\rho$ measures the number of decimal digits of accuracy gained in one iteration (for sufficiently large $N$)\n", "\n", "
Example (cont.). \n", "\n", "Let's go back to our example: $f(x) = e^x - 2x - 1$. We have defined $g(x) = \\log( 2x + 1 )$ and we are trying to solve the fixed point problem $g(x) = x$ (which is mathematically equivalent to $f(x) = 0$).\n", "\n", "We have seen (numerically) that $x_{n+1} = g(x_n)$ with $x_1>0$ converges to $\\xi = g(\\xi) \\approx 1.26$. What is the order of convergence? What is the asymptotic error constant? Asymptotic rate of convergence?\n", "
" ] }, { "cell_type": "code", "execution_count": 11, "id": "d01a6a42", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "μ (generic function with 2 methods)" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function μ( x, ξ, α=1 )\n", " return @. abs( x[2:end] - ξ ) / ( abs(x[1:end-1] - ξ )^α );\n", "end " ] }, { "cell_type": "code", "execution_count": 12, "id": "736d24e9", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "21-element Vector{Float64}:\n", " 0.999205051893074\n", " 0.99841359441868\n", " 0.996841052928944\n", " 0.9937367881169422\n", " 0.9876863189309779\n", " 0.97617878840982\n", " 0.9552662983832602\n", " 0.9201668651834529\n", " 0.8682070800166984\n", " 0.8034158099852159\n", " 0.7370073130661621\n", " 0.6804438941331454\n", " 0.6388424557889528\n", " 0.6111968825113593\n", " 0.5939641345990201\n", " 0.5836258577310153\n", " 0.5775606878797853\n", " 0.5740478904522482\n", " 0.5720282875181989\n", " 0.5708720331145754\n", " 0.5702116439591833" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# this is e_{n+1}/e_{n} for n=1,...,14\n", "μn = μ( x, ξ)\n", "ρ = -log10( μn[end] );\n", "\n", "μn" ] }, { "cell_type": "markdown", "id": "08dfca60", "metadata": {}, "source": [ "
\n", "\n", "This quantity appears to converge to $\\mu \\in (0,1)$. Therefore, we have $(x_n) \\to \\xi$ linearly with asymptotic error constant $\\mu$ and asymptotic rate of convergence $\\rho$ where these constants are as follows:\n", "\n", "
" ] }, { "cell_type": "code", "execution_count": 13, "id": "b91af37e", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(ξ, μn[end], ρ) = (1.2564312086261678, 0.5702116439591833, 0.24396391846186502)\n" ] } ], "source": [ "@show (ξ, μn[end], ρ);" ] }, { "cell_type": "markdown", "id": "9cbebe89", "metadata": {}, "source": [ "
\n", "\n", "Next: we will now prove (mathematically) that $x_{n+1} = g(x_n)$ converges linearly with $\\mu = |g'(\\xi)|$. Our computed value of $\\mu$ (after $15$ iterations) is close to this theoretical value: \n", "\n", "
" ] }, { "cell_type": "code", "execution_count": 14, "id": "6daf0600", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "g_prime(ξ) = 0.5693362740816775\n", "μn[end] = 0.5702116439591833\n" ] }, { "data": { "text/plain": [ "0.000875369877505805" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@show g_prime( ξ )\n", "@show μn[end]; \n", "abs( g_prime( ξ ) - μn[end] )" ] }, { "cell_type": "markdown", "id": "23b03b7d", "metadata": {}, "source": [ "* Contraction: $g : [a,b] \\to \\mathbb R$ is a contraction if there exists $L \\in (0,1)$ such that $|g(x)-g(y)|\\leq L|x-y|$ for all $x,y \\in [a,b]$," ] }, { "cell_type": "markdown", "id": "ee6a08d2", "metadata": {}, "source": [ "
Contraction Mapping Theorem (Banach fixed point theorem). \n", "\n", "Suppose $g:[a,b] \\to[a,b]$ is a contraction. Then, there exists a unique fixed point $\\xi = g(\\xi) \\in [a,b]$. Moreover, the iteration $x_{n+1} = g(x_n)$ converges to $\\xi$ for all $x_1 \\in [a,b]$\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "d8b91a40", "metadata": {}, "source": [ "
Proof. \n", "\n", "By Brouwer's fixed point theorem ($g:[a,b] \\to [a,b]$ is continuous), there exists a fixed point $\\xi \\in [a,b]$. If $\\zeta$ is another fixed point, then there exists $L \\in (0,1)$ such that\n", "\n", "\\begin{align}\n", " |\\xi - \\zeta| &= |g(\\xi) - g(\\zeta)| \n", " \\tag{as $\\zeta = g(\\zeta)$}\\nonumber\\\\ \n", " %\n", " &\\leq L |\\xi - \\zeta| \n", " \\tag{as $g$ is a contraction} \n", "\\end{align}\n", "\n", "That is, $0 \\leq (1-L) |\\xi - \\zeta| \\leq 0$ and so $\\xi = \\zeta$. i.e. $\\xi$ is unique.\n", "\n", "Defining $x_{n+1} = g(x_n)$ for $x_1 \\in [a,b]$, we have \n", "\n", "\\begin{align}\n", " |x_{n+1} - \\xi| &= |g(x_n) - g(\\xi)| \\nonumber\\\\\n", " %\n", " &\\leq L |x_{n} - \\xi| \\nonumber\\\\\n", " %\n", " &\\leq L^n |x_1 - \\xi|.\n", "\\end{align}\n", "\n", "Since $L \\in (0,1)$, we have $L^n \\to 0$ as $n \\to \\infty$ and thus $(x_n) \\to \\xi$.\n", "
" ] }, { "cell_type": "markdown", "id": "8a427429", "metadata": {}, "source": [ "
Remark. \n", "\n", "Notice that, by the mean value theorem, there exists $\\eta_n$ between $x_{n}$ and $\\xi$ for which \n", "\n", "\\begin{align}\n", " \\frac\n", " {|x_{n+1} - \\xi|}\n", " {|x_n - \\xi|}\n", " %\n", " &= \\frac\n", " {|g(x_n) - g(\\xi)|}\n", " {|x_n - \\xi|} \\nonumber\\\\\n", " %\n", " &= \\frac\n", " {|g'(\\eta_n)| |x_n - \\xi|}\n", " {|x_n - \\xi|} \\nonumber\\\\\n", " %\n", " &= |g'(\\eta_n)|.\n", "\\end{align}\n", "\n", "Therefore, since $(x_n) \\to \\xi$, and $\\eta_n$ is between $x_n$ and $\\xi$, we have \n", "\n", "\\begin{align}\n", " \\lim_{n\\to\\infty}\n", " \\frac\n", " {|x_{n+1} - \\xi|}\n", " {|x_n - \\xi|}\n", " %\n", " &= |g'(\\xi)|\n", "\\end{align}\n", "\n", "and the order of convergence is at least linear with asymptotic error constant $\\mu = |g'(\\xi)|$. \n", "\n", "(Recal that, if $\\mu = 0$, the order of convergence is superlinear)\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "11fd62e9", "metadata": {}, "source": [ "
Example (cont.). \n", "\n", "Show that $g(x) = \\log( 2x + 1 )$ is a contraction on $[1,2]$ with $L = 2/3$. \n", "Therefore $x_{n+1} = g(x_n)$ converges for all $x_1 \\in [1,2]$.\n", "\n", "\n", "Answer: $g' \\in [\\frac23, \\frac25]$ on $[1,2]$ and so $g$ is a contraction with constant $L = \\frac23$. Applying the contraction mapping theorem shows that $(x_n) \\to \\xi$ in this case.\n", "
" ] }, { "cell_type": "markdown", "id": "b891928d", "metadata": {}, "source": [ "* How many iterations to guarantee an error of $ \\leq \\epsilon$?\n", "\n", "Recall that $|x_{n+1} - \\xi| \\leq L |x_{n} - \\xi|$. In particular, we have $|x_1 - \\xi| \\leq |x_1 - x_2| + |x_2 - \\xi| \\leq |x_1 - x_2| + L|x_1-\\xi|$ and thus $|x_1 - \\xi| \\leq \\frac1{1-L} |x_1 - x_2|$. As a result, we have \n", "\n", "\\begin{align}\n", " |x_{n+1} - \\xi| &\\leq L^n |x_1 - \\xi| \\nonumber\\\\\n", " &\\leq \n", " \\frac\n", " {L^n}\n", " {1-L} \n", " |x_1 - x_2|\n", "\\end{align}\n", "\n", "This number is less than $\\epsilon$ when\n", "\n", "\\begin{align}\n", " n \\geq \\frac{ \\log \\frac{|x_1-x_2|}{\\epsilon(1-L)} }{\\log\\frac1L}.\n", "\\end{align}\n", "\n", "This gives an upper bound on the number of iterations needed in order to be within some fixed tolerance.\n", "\n", "
Example. \n", "\n", "Back to $f(x) = e^x - 2x - 1$ written as the fixed point problem $g(x) = x$. For an error of $\\epsilon = 10^{-10}$, we require around $54$ iterations. Actually, $40$ will do (but we have an upper bound): \n", "
" ] }, { "cell_type": "code", "execution_count": 15, "id": "bb7d1788", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "N = log(abs(x1 - x2) / (ϵ * (1 - L))) / log(1 / L) = 53.78490871074893\n" ] }, { "data": { "text/plain": [ "8.923284333661741e-11" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ϵ = 1e-10; \n", "L = 2/3;\n", "x1 = 1; x2 = g(x1);\n", "@show N = log( abs(x1-x2) / ( ϵ * (1-L) ) ) / log( 1/L )\n", "\n", "\n", "x = simple_iteration( g, 1.0, 54, 1e-20)\n", "\n", "abs( x[40] - ξ )" ] }, { "cell_type": "markdown", "id": "a73df67f", "metadata": {}, "source": [ "
Definition. \n", "\n", "$\\xi$ is a *stable fixed point of $g$* if the iteration $x_{n+1} = g(x_n)$ converges to $\\xi$ for all $x_1$ sufficiently close to $\\xi$ \n", "$\\xi$ is an *unstable fixed point of $h$* if for all $\\delta > 0$ there exists $x_1 \\in [\\xi-\\delta, \\xi+\\delta]$ for which $x_{n+1} = g(x_n)$ does not converge to $\\xi$.\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "f16a9887", "metadata": {}, "source": [ "
Example (cont.). \n", "\n", "$\\xi = 1.26...$ is a stable fixed point of $g(x) = \\log( 2x + 1 )$. \n", "Proof: contraction mapping theorem on $[1,2]$\n", "
" ] }, { "cell_type": "markdown", "id": "dd725c3f", "metadata": {}, "source": [ "* Theorem: Suppose $h(\\xi) = \\xi$ and $h'$ is continuous with $|h'| > L > 1$ near $\\xi$. Then, $x_{n+1} = h(x_n)$ does not converge to $\\xi$ for all $x_1 \\not= \\xi$. " ] }, { "cell_type": "markdown", "id": "52c2225b", "metadata": {}, "source": [ "
Proof. \n", "\n", "Let $I_\\delta = [\\xi-\\delta,\\xi+\\delta]$ be such that $|g'| > L$ on $I_{\\delta}$. If $x_n \\in I_{\\delta}$ then there exists $\\eta \\in I_\\delta$ such that\n", "\n", "\\begin{align}\n", " |x_{n+1} - \\xi| &= |h(x_{n}) - h(\\xi)| \\nonumber\\\\\n", " %\n", " &= |h'(\\eta)| |x_n - \\xi| \\nonumber\\\\\n", " %\n", " &> L |x_{n} - \\xi|.\n", "\\end{align}\n", "\n", "Similarly, if $x_{n+1} \\in I_\\delta$, then $|x_{n+2} - \\xi| > L^2|x_{n} - \\xi|$. Since, $L> 1$ (and so $L^n\\to\\infty$ and $n\\to\\infty$), there exists $N$ such that $x_{n+N} \\not\\in I_{\\delta}$. \n", "\n", "We have therefore seen that if $x_n \\in I_\\delta$, there exists $N$ such that $x_{n+N} \\not\\in I_{\\delta}$. As a result, $(x_n)$ does not converge to $\\xi$. \n", "\n", "
" ] }, { "cell_type": "markdown", "id": "1c3ca153", "metadata": {}, "source": [ "
Example (cont.). \n", "\n", "$\\xi = 1.26...$ is an unstable fixed point of $h(x) = \\frac{e^x -1}{2}$. In fact, for $x_1 > \\xi$, $x_{n+1}= h(x_n)$ diverges to infinity. \n", "$\\zeta = 0$ is a stable fixed point of $h$. For $x_1 < \\log 2$, $x_{n+1}= h(x_n)$ converges to $0$. \n", "
" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.7.1", "language": "julia", "name": "julia-1.7" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.7.1" } }, "nbformat": 4, "nbformat_minor": 5 }